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ABSTRACT 

The influence of optical scattering and thermal radiation models an the Yarkovsky-O'Keefe- 
Radzievskii-Paddack (YORP) effect is studied. Lambertian formulation is compared with 
Hapke scattering and emission laws and Lommel-Seeliger reflection. Although the form of 
reflectivity function strongly influences mean torques due to scattering or thermal radiation 
alone, their combined contribution to the rotation period YORP is not much different from the 
standard Lambertian values. For higher albedo values the differences between the Hapke and 
Lambert models become significant for the YORP in attitude. 
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1 INTRODUCTION 

The importance of radiation recoil forces on both orbital motion 
and rotation of minor bodies in the Solar System has been com- 
monly appreciated over the last decade. The Yarkovsky effect, 
caused by lagged thermal radiation from the surface of a spin- 
ning body (directly d etected in the orbital motion of 64 89 Golevka 
dCheslev et ai1l2003T) and 1992 BF dVokrouhlickv et alj|2008h ), has 
occurred to be the key to the proper understanding of asteroid 
long-term dynamics. Since the paper of iRubincaml ( 120000 . the in- 
fluence of torques due to radiation recoil is known as the YORP 
( Yarkovsky-0 ' K eefe-Radz i evskii- P addack) effec t , ackn o wledging 
the w orks of lYarkovskvl dl90lh . iRadzievskiil dl954T) , IPaddackl 
dl969h . and lO'Keefeldl976h . Unlike the orbital Yarkovsky effect, 
YORP involves both the scattering of incident light and its thermal 
re-radiation, and it occurs even for objects with zero conductiv- 
ity. Direct det ections of YORP effect in the rotatio n of asteroids 
54509 YORP dLowrv et alj200l iTavlor et al. | l2007I). 1862 Apoll o 
dKaasalainen et alj|20071). 1620 Geo graphos dDurech et ai] ;2008). 
and 3103 Eger dDurech et al.l 120090 . has proved the existence of 
YORP. On the other hand, however, the agreement between the ob- 
served and modeled values in each case can be qualified as merely 
having a similar order of magnitude and all present YORP mod- 
els are still simplified and incomplete. What is worse, the failure 
to detect a theoret ically predicted YORP effect for 25143 Itokawa 
dDurech et alll20080 has helped to realize an extreme sensitivity of 
these simplified models to the fine details of shape, c entre of mass 
location and spin axis orientation in the body fram e (Durech et al. 
l2008tlScheeres & Gaskellll2008l ; IStatlenl200alBreiter et alj2 009). 

Most of these models assume that both scattering and ther- 
mal radiation is Lambertian, i.e. a photon can be emitted or scat- 
tered with equal probability in any direction, hence the exiting 
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flux is proportional to the cosine of zenith distance according to 
the pr ojected area of the ra diating/scattering surface element. Al- 
though iBreiter et alj ([2007) mentioned a more general scattering 
model, their work on YORP for spheroids was based on the Lam- 
bertian assumption. IScheeresI (2007) added a specular reflection to 
the Lambertian scattering, but such an improvement h as no effect 
on the mean torque, as obser ved by Rubincam d2000Q and proved 
by iNesvornv & V okrouhlickv (2008) and Rubi ncam & Pa ddack 



d2010P . IStatleir ?2009) made a step further, using a hem ispheric 



albedo derived from the scattering model of Hapke (2002) instead 
of the usual Bond albedo ( I Vokrouhlickv & Bottkdl200lh . More- 
over, the TACO model of Statler for the first time incorporates the 
important observation that photons bouncing between various sur- 
face elements do not produce the net torque until they they finally 
exit into outer space. 

The main objective of the present work is to include a non- 
Lambertian scatte ring and radiation into the recent YORP model of 
IBreiter et al J J201oh and judge the significance of this improvement. 
Roughly speaking, a departure from the Lambertian model is essen- 
tially caused by interreflections and occlusions. Both phenomena 
occur at various levels of resolution and one has to be careful about 
this issue. Out of the two principal scattering models for a steroids 
surface developed by Lumme & Bowell ( 19 8l|) and Hapke dHapkd 
1 19811 : lHapke & Wellslll98ll ; lHapkelll984 1 19861 12003. l2008h we 
have chosen the latter. However, both models were created to in- 
terpret photometric observations; as such, they attempt to include 
phenomena happening at various resolution levels that merge in 
the final integrated brightness. In these circumstances, our present 
paper focuses on accounting for the regolith grain size (< 1 mm) 
scale phenomena described by an appropriate part of the Hapke re- 
flectance and emissivity models. It means ignoring the macroscopic 
roughness corrections and shape-dependent beaming factors. In- 
terreflections occurring between larger surface fragments require a 
different approach and will be discussed in another article, whereas 
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the large scale occlusions (shadowing) are already incorporated in 
most of existing YORP models. To a large extent the present con- 
tribution has been motivated by the problem of YORP effect on a 
high albedo asteroid 3103 Eger, where only a convex shape model 
is available, so larger scale interreflections have to be neglected 
anyway. 

We decided to present a detailed derivation of the ra- 
diation recoil force and associated YORP torque using the 
terminology of modern radiometry instead of traditional 
astrop hysical framework dating back to Chandrasekhar 
( 1950). In this respect we owe much to the collection of 
Max Fairbairn essays available on-line thanks to J.B. Tatum 
(http : // orca . phys . uvic . ca/" tatum/plphot . html). 



2 SCATTERING AND RADIATION 
2.1 Irradiance 

Consider an infinitesimal element dS of an atmosphereless celestial 
body surface. Further, we call dS a physical surface, to distinguish 
it from a normal surface i.e. an infinitesimal surface perpendicular 
to some specified direction of incident or emitted radiation. 

Local solar frame (LSF) will be defined with the origin at the 
centre of dS , with x-axis pointing to the intersection of the merid- 
ian passing through the Sun and the horizon plane, z-axis directed 
along the outward normal to the physical surface (i.e. to the local 
zenith), and y-axis completing the right-handed orthogonal system. 
Then the unit vector directed to the Sun has a simple form 



where fi is the cosine of the zenith distance and (p is the azimuth 
angle in LSF. The cosine of the phase angle between s and o will 
be designated //, and defined as the scalar product 



A' 



■ O = Sq yj 'l -fi 2 COS (p + p Q jJ. 



(6) 



The power flux d 2 <I> r scattered from dS into a solid angle dQ. 
in direction o is described by reflected radiance L r 



d 2 O r 



(7) 



pdSdQ 

where pdS is the normal surface perpendicular to o. Writing L r (o) 
we should bear in mind an implicit dependence on the Sun direction 
s, because the reflected power depends on the incident flux from 
the Sun as well. This dependence becomes more explicit, when we 
introduce a bidirectional reflectance distribution function (BRDF) 
f r defined as the ratio of the radiance L r reflected in the direction o 
to the irradiance E from the energy source located in the direction 
s 



Mi, o) ■■ 



Lr(Q) 
E(s)- 



Although a bidirectional reflectance (BDR) function p 
p(s,o)=fi @ Ms,o), 



(8) 



(9) 



seems to be more common in planetary photometry than BRDF, we 
choose f t as a more convenient quantity offering, for example, the 
reciprocity relation / r (s, o) = f(o, s). Using Eq. {8]l we can express 
the reflected radiance as 



So 


Ho 



(1) 



depending only on the solar zenith distance through its cosine /j 
and sine 



so ■■ 



(2) 



If the Sun is located at the distance R & , the collimated radia- 
tion flux density (power per normal area perpendicular to s) arriving 
from the point R & s is 

\2 



J = Jo 



(3) 



where the Solar constant Jq ss 1366 WirT 2 is defined for nominal 
distance Ro = I au. Irradiance or insolation E of an arbitrarily ori- 
ented surface element is the ratio of the incident power flux dO; 
to the physical area dS . Accounting for the area projection factor 
s-n = Ho, where n is the unit vector directed to zenith, we can write 

dOj 
~dS 



E(S): 



vjp. . 



(4) 



The visibility function v is either 1, when the Sun is visible at dS, 
or 0, when the Sun is occluded. 



2.2 Bidirectional scattering 

Incident radiant power dO; = EdS is partially absorbed (converted 
into heat) and partially reflected in various directions o 



yU 2 COS(j 

•yu 2 sin if 



(5) 



^r(0) = f r (s,0)E(s) = vfl-(s,0)fl Q J. 



(10) 



Recalling for the reference a traditional, Lambertian BRDF 
with albedo A 



,/L — ~~> 
n 



(11) 



we adopt the anisotropic BRDF proposed by Hapke, namely its 
version from lHapkel l2002l) without macroscopic roughness effects 



■[(1+B)P + H(ji q )H(m)-1], 



f(s,o)= 

47r(p +yu) 

where the Henyey-Greenstein particle phase function is 

/> = (l-g 2 )[l + 2 g/ / + /]^, 

and the opposition surge function B is defined as 

-1 



S = S 
with 



h V 1 -ju' 



fin 



w (l-g) 



(12) 



(13) 



(14) 



(15) 



The Chandrasekhar multipl e scattering function H is defined 
in terms of an integral equation dChandrasekhari ri950), but we use 
the explicit second order approximation given bv lHapkel ( 120021) 



H(x) = 
where 



, , l-2r x l+x 
1 - wx \r + In 

2 x 



(16) 
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l- VT^ 

1 + vr^; 



(17) 



Thus, apart from the incoming and scattered flux directions s and 
o, the Hapke BRDF depends on four physical parameters of the 
surface: single scattering albedo w, regolith compaction parameter 
h, the opposition surge amplitude Sq (sometimes replaced by So), 
and the asymmetry factor of the Henyey-Greenstei n function g. A 
more recent version of reflectance was proposed bv lHapkel l2008); 
the modification amounts to adding the dependence on porosity as 
a multiplicative factor in / r and a divisor in the argument of H. 
This modification is easy to implement, but we suspend its use 
until controversies concerning the depende nce of the opposition 
effect on porosity are resolved (Hapke 2008). 

The total power flux O e emitted into the hemisphere 
fi+ = [(jd,(f>) : O«ju«1,O«0<2ct}, 
divided by the emitting physical area is termed radiant exitance M 
d* e 



M ■ 



dS 



(18) 



Recalling the definition Q, we find for the exitance due to scatter- 



(19) 



f d 2 o r r 

J n+ dS dQ, J n+ 

But, according to Eq. l[8}, radiance is related to irradiance by the 
BRDF f r , hence 



MM = E(s) f / I (s,o) J udQ, 
Jn + 



(20) 



and the dependence on the Sun location s appears explicitly. 

At this point we can introduce the notion of hemispheric 
albedo A/, as the ratio 

My 



E(sY 



(21) 



Combining Eqs. l !20t and J2U we see that 

A,0"o) -■ 



/ r (s, o)ndCl= d(p f t (s,o)ndn, (22) 
n+ Jo Jo 



and Eq. i20\ is simplified to 

M t (s) = E{s)A h (p Q ). (23) 

For a given set of Hapke parameters the integral l !22t can be 
evaluated numerically on a sufficiently dense set of p Q values, al- 
lowing to construct an appropriate approximating function. Using 
the least squares adjustment we construct 

1 - «2y"o 



p & A h (p ) ss A B p + ain& - 

1 + ayiQ 

where the Bond albedo Ab, defined as 



= - f A h (p Q )p Q dn Q = 2 \ A h (p Q )p Q dp Q , 
" Jn + Jo 



(24) 



(25) 



is the mean slope of the product yu Ah, and the coefficients a,- of a 
simple rational approximation describe the deviation from the lin- 
ear model. We focus on the properties of yU Ah, because in next 
sections the hemispheric albedo will always appear multiplied by 
the cosine of the Sun's zenith distance. Note, that the adjustment 
of A^(jj q ) leads to different values of a,-, degrading the quality of 
approximation of the product /J© Ah. 



2.3 Geometric albedo 

Although the geometric (or physical) albedo is not directly involved 
in the computation of radiation recoil force, we need it to select 
an appropriate value of w for the Hapke model, since usually the 
observations provide for an asteroid only the geometric albedo and 
the spectral type. 

Let us begin with the notion of intensity I. In contrast to the 
previously discussed quantities, intensity refers to the power dO 
emitted from the whole body surface (not only from an infinitesimal 
dS) in some direction q, divided by the solid angle dCl centered at 
9 

. d<D 

/(?) = ^r. (26) 
dQ 

The geometric albedo p is the ratio of observed intensity of some 
presumably spherical object to the intensity of Lambertian disk 
with the same diameter as the assumed sphere - both observed from 
the direction to the Sun, i.e. with zero phase angle. This leads to the 
integral definition 



p = 2n p Q f r (s,s)dp B . 



(27) 



IVerbiscer & Veverki] dl995h provide expressions that allow to 
compute Hapke parameters h, Bq, g for various spectral types as 
functions of a given geometric albedo p and the mean slope param- 
eter G of the IAU two-parameter magnitude system jBowell et al.l 
Il989h . 

2.4 Directional thermal emission 

The energy leaving a surface element dS does not consist only of 
scattered radiation. If the element has temperature T > 0, it also 
emits thermal radiation. Radiant exitance M\, through Q.+ for a 
black body is given by the Stefan-Boltzmann law 

dO b 



dS 



■<tT\ 



(28) 



where ar = 5.67 x 10~ 8 WirT 2 KT 4 , is the Stefan-Boltzmann con- 
stant, and Ob is the black body value of a more general thermal 
radiation power flux O t - The point black body radiation is by defi- 
nition isotropic, whereas a black body surface radiation is Lamber- 
tian, so the associated radiance Lb(o) in the direction o is obtained 
from the general definition of a thermally emitted radiance L t anal- 
ogous to < fT0b 

d 2 O t 



(29) 



(idSdQ, 

dividing the exitance Mb by the '^-averaged' solid angle of a hemi- 
sphere n 

L b = ^. (30) 

71 

Indeed, using J30I ) and the definition of exitance, we verify that 

f L b //dn = M b . (31) 

Jn + 

Hemispheric emissivity eh is defined as the ratio of actual ther- 
mal exitance Af t to the black body exitance Mb 

M, 

* = w b (32) 

This global quantity should not be confused with a directional 
emissivity e(o), defined as the ratio of radiances 
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e(o)- 



U(o) 
L h (o) 



(33) 



Directional emissivity plays the role similar to that of BRDF in 
scattering, although their definitions essentially differ: the former is 
dimensionless ratio of two radiances, while the latter (with dimen- 
sion sr _1 ) is the ratio of radiance to irradiance. So, the thermally 
emitted radiance in the direction of o can be expressed as 



L t (o) = e(o)L b (o)=^-o-T 4 . 

71 



(34) 



In pr esent study we adopt the dir ectional emissivity function 
of Hapke dHapkefl 993; Lag errosll 19961) 



= Vl-w/%), 



(35) 



where w is Hapke's single scattering albedo, and the Chandrasekhar 
function H is given by Eq. fll6t . Thus the emitted radiance is 



L t (o) = e(ji) 



crT 4 Vl-w „, , „ 4 



■H( f i)crT\ 



(36) 



7X n 
The exitance M t resulting from Eq. d36t is 

M t = I Lt(o)ijdfl. 
Jn + 

Confronting it with the primary definition of d32t . we can use 
the relation 



(37) 



M t = e h crT* 



(38) 



which leads to the integral expression of hemispheric emissivity 

(39) 



Jn + n Jo 



evaluating to a single number for a given set of Hapke parameters. 

The infrared radiation of asteroids is often related with the no- 
tion of beaming effect, empirically accounted for by the beaming 
factor r\ dLebofskv & S pencer 1989; Lagerros 1993). We do not in- 
troduce the beaming factor in our model for a number of reasons: 
i) the part of beaming that depends on grain size scale radiation 
transfer should be present in the emissivity function of Hapke, ii) 
the contribution of thermal lag to the beaming factor is present in 
the surface temperature model with conductivity , iii) larger scale 
radiation exchange contribution l lLagerrosIl 19981) will be included 
in future extensions of our model together with optical interreflec- 
tions. 



2.5 Energy balance 

Conservation of energy implies, that the total power scattered, ther- 
mally re-emitted, and conducted inside the body, should be equal to 
the incident power flux <!>;. In terms of power density (per physical 
surface) it means that 



E(s) = M I (s) + M t -Q, 



(40) 



i.e. irradiance E is equal to the sum of total radiant exitance M and 
of the conducted heat flux density (-0. Given a nonzero surface 
conductivity K, we have 



Q = -KnVT, 
and then 

E(s) = A h Ou ) E(s ) + 6 h o- T 4 + K n ■ VT, 



(41) 



(42) 



e h o-r 4 = v^ (l-A h Ou ))/+<2. 



(43) 



If K = 0, Eq. J43t directly provides the surface temperature, gen- 
eralizing the usual Lambertian Rubincam approximation of the 
YORP effect to the Hapke refiectance/emissivity model. With K + 
0, Eq. J43t serves as a boundary condition for the heat conduction 
problem. 



3 RADIATION RECOIL FORCE AND TORQUE 
3.1 Force expression 

Photon flux, leaving dS in direction o, carries energy and momen- 
tum (energy divided by the velocity of light c), inducing the recoil 
force F equal to the time derivative of momentum and directed op- 
posite to o. The force can be easily expressed in terms of emitted 
radiance, provided we introduce a radiance vector 



L(o) = L(o)o, 



(44) 



where L is the sum of scattered and thermal L t . The definition of 
radiance implies that the force density in the direction o per physi- 
cal area and solid angle is 

d 2 F d 2 (* r + cD t )o 



-=-^L(o). 

c c 



(45) 

dSdQ. dSdQ c c 

Integrating over the hemisphere Q.+, we find the net force per phys 
ical area 



■ 7F =— f fiL(o)dCl. 

c J n+ 

Substituting the expressions d 1 Ob and ( 134b . we have 

or, observing the independence of directional emissivity on azimuth 

2o-T 4 n rl 



(46) 



(47) 



dF _ E(s) C 



fif I (s,o)odfl- 



Jo 



^c(fi)dfi. 



(48) 



Using LSF we can conveniently decompose the force density 
into the sum of two perpendicular components along the axes z and 
x, i.e. along the surface normal n and the unit vector m 



(49) 





' ' 




( 1 ' 


n = 





, m = 







i 1 , 




. o J 



(50) 



In an arbitrary reference frame we can compute m as 

»i = (*-^o")^o 1 , 

taking //© = sn. 

Splitting the first integrand in Eq. J48 1 > into a sum 

/j f r (s,o)o = frfj 2 n + fail yjl - jj 2 cos(pm 

+ftH yj 1 -yU 2 sin</> (nxm), 

and recalling that f r (s,o) is an even function of azimuth <p, hence 
the last term in Eq. (51) is odd and its integral over Q+ does vanish, 
we introduce two auxiliary functions 

h(Ho) = | d<f> I ^Ms^dfi, (52) 



(51) 



f d<j> f fl Z f r (s,0)dfl, 

Jo Jo 

f d(f> f fiJl-fj 2 cos(pf [ (s,o)dfi, 

Jo Jo ' 



(53) 



Table 1. Sample Hapke parameters and related quantities 





p = 0.6, type hj 


p — 0.1, type S 


w 


0.856 


0.139 


h 


0.044 


0.049 




0.8576 


1.5407 


S 






Ab 


0.47542 


0.046712 


ai 


0.14550 


0.031021 


"2 


1.5880 


1.6968 


a 3 


0.9741 


3.6800 


fio 


0.09590 


0.01212 


fit 


-0.03374 


0.3383 




2.1268 


1.4786 


fill 


0.032156 


0.006259 


61 


0.42458 


0.16050 




0.03024 


0.00255 




0.58832 


0.96905 



as well as a coefficient 

h = 2 f //ef^d//, 
Jo 

that allow to rewrite the formula ((48} as 

dF E(s) r n / 3 crr 4 

— = [IiQi )n + I 2 (Me)m\ n. 

dS c c 

Substituting the boundary conditions J42t into d55t we remove 
the explicit dependence on T 4 , obtaining 



(54) 



(55) 



dF 

dS 



m 

c 



j , , , , l-AhQi ) 1 
eh 



6/3 

ceh 



n. (56) 



However, we prefer to rearrange the force expression into a more 
comprehensive form 



dF 2 1+f . vJ 

— = -- -[vJpi Q + Q\n+ — [Xin-X 2 m], 

dS 3 c c 



(57) 



where functions arguments have been omitted for the sake of 
brevity. 

The coefficient 

3 h 



2 6 h 

is a small quantity of order 1CT 2 or less. The two functions 
A h I 3 



ft 



—h 



f*eh, 



(58) 

(59) 
(60) 



also represent a small deviation from the Lambertian model (see 
Fig.Q]based upon the data from Tab.[T|l. 

In our judgement there is no point in producing excessively 
accurate approximations of corrections to the Lambertian model, 
so we use relatively simple functions, found by trial and error, 



Xi ~ %WH<i 



X 2 ~ flOAD (3-^21^0+^q) Jj7 



1-/43 



(61) 
(62) 



with coefficients generated by the least squares adjustment to the 
results of numerical quadratures. 
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type E, p=0.6 




type S, p=0.1 



0.004 - 



0.002 



-0.002 - 




Figure 1. Functions describing the non-Lambertian force model. Solid line 
represents the remainder PeiA^ - Ab), dashed and dotted lines refer to X\ 
and X2 respectively. 



The limit case of Lambertian model results from setting f = 
X\ = X2 = 0, and then Eq. d57| l simplifies to 



dF l 
dS 



3c 



(vJ(i Q + Q)n. 



(63) 



Of course, this step also requires assuming a constant = Ag in 
boundary conditions J42b for a heat conduction solver providing Q. 



3.2 Torque expression 

The force defined by Eq. \57\ generates for each surface element a 
torqu^] 

/ dF\ 2 \+£ T 

AM = rx— dS =--_*.[v/ P0 + 0](rxdS) 



+— [Xi (rxdS)-X 2 (rxmdS)]. 

c 



(64) 



The two cross products present in this formula differ in nature; the 
first, namely 

rxdS = rxndS, (65) 
is constant over time in the body frame, whereas the second, 

1 We maintain the symbol M from previous papers hoping that it will not 
be confused with exitance M appearing only in Sect. [2] 
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rx s - iir,rxn rxs u n 
rxmdS = — AS = AS- — (rxdS), (66) 

% Sq Sq 

is time dependent due to the solar motion on local celestial sphere 
ofdS. 



3.3 YORP effect computation 

Total YORP torque M resulting from the force d57b 



Table 2. Physical and orbital data for the test bodies 



M 



AS, 



(67) 



is obtained by integration over the body surface. The way the in- 
tegration is handled depends on the type of a body shape model: 
it can be performed analytically if the surface equation is explicit 
(e.g. spherical harmonics expansion) or, more often, replaced by 
the sum over flat faces of a triangulation mesh. In the Rubincam 
approximation, when 2 = 0, one may simply substitute Eq. l !57t 
into l !67t to obtain the torque for a given position of the Sun in the 
body frame. Most often the resulting torque values are then aver- 
aged with respect to rotation and orbital motion in order to extract 
the secular effects in rotation rate and attitude dynamics. This step 
requires assumptions about the nominal rotation model that pro- 
vides the averaging kernel and solar ephemerides. 

When the heat conduction is included, the nominal rotation 
model enters much earlier than in the final averaging: surface 
temperature oscillations are lagged with respect to the insolation, 
hence we cannot find Q, required by the torque formula, without 
the knowledge of rotation history. Choosing the simplest princi- 
pal axis rotation mode (known as the gyroscopic approximation) 
we can easily add non-Lambertian corrections to the algorithm of 
IBreiter et alJd2010h based on a nonlinear ID thermal model. 

The one-dimensional model, where conduction is restricted to 
the direction normal to the surface, allows a separate treatment of 
each triangular face of the shape mesh. There, having specified the 
obliquity e (angle between the spin axis and normal to the orbit) 
we sample mean anomaly and rotation phase creating the vector of 
absorbed radiant flux values (the first term in the right-hand side of 
Eq. (143b V Its discrete Fourier transform (DFT) serves to compute 
the DFT spectrum of Q by an iterative process. Once the spectrum 
of Q is known, one is able to compute the torque M. Until this step, 
no essential modifications of the algorithm are required; all one has 
to do is to replace the constant albedo A ( understood as the B ond 
albedo Ag) in the boundary conditions of IBreiter et alj ( 120101) by 
the hemispheric albedo function Ah(p e ). Apart from this point, the 
heat diffusion solver remains practically unchanged; however, com- 
puting the mean torque demands a deeper revision. The spectrum of 
absorbed power flux v.//j©(1 - Ah), evaluated for the conductivity 
contribution, cannot be recycled in the Rubincam part, because the 
hemispheric albedo is not a constant. Thus, outside the conductivity 
related block, we directly compute the mean values of projections 
of the Rubincam part AM (given by Eq. d64t with Q = 0) on unit 
vectors 



e\ = sinfi'e A +cos£2'e v , 
e2 = - cosCl' e x + sinCl' e v , 
e 3 = e z> 



(68) 



where e x , e y , e z form the body-fixed frame basis, a nd Q! is the 
rotati on phase measured from the asteroid's equinox ( Brei ter et al.1 
prime added to avoid confusion with solid angle Q. of the 
present paper). Then we add the mean values resulting from the 
DFT spectrum of Q, obtaining the final averaged torque projections 
{Mi) = {M- ei ). 







(54509) YORP 


(3103) Eger 


epoch 


JD 


2452117.5 


2446617.0 


semi-axis 


au 


0.9930 


1.4068 


eccentricity 




0.2305 


0.3548 


inclination 


deg 


1.9971 


20.939 


asc. node 


deg 


283.835 


129.972 


arg. perihelion 


deg 


272.091 


253.661 


rotation period 


h 


0.2029 


5.7102 


ecliptic pole (A,/}) 


deg 


(180,-85) 


(224, -72) 


effective diameter 


m 


113 


1778 


density 


kgirT 3 


2500 


2800 


conductivity 


Wm-'K" 1 


0.02 


0.02 


specific heat 


Jkg-'K" 1 


680 


680 


max. mom. inertia 


kgm 2 


3.04 xlO 12 


2.80 x 10 20 



If u> stands for the rotation rate, the dynamics in the gyroscopic 
approximation is governed by 

Mi 

E = — F* 

M 2 



(69) 
(70) 
(71) 



coC tans 

where C designates the maximum moment of inertia in the princi- 
pal axes frame. 

The conclusion of IBreiter et al.1 ( l2010h . that all ID thermal 
models imply the independence of the mean period related com- 
ponent (M3} on conductivity, holds true regardless of the scattering 
and emission laws. 



4 EXEMPLARY RESULTS 

In order to see how the improvement of scattering and emission 
laws affects the simulated YORP effect, we considered two exem- 
plary objects out of the four known to have observationally con- 
firmed spin acceleration: (54509) YORP asteroid wit h an irregu- 
lar, radar-determined shape modePI ( Taylor et al J2007I) . and (3103) 
Eger with a convex shape model obtained by lightcure inversion 
jDurech et al.l2009f) . Both shape models, consisting of 572 (YORP) 
and 1972 (Eger) triangular facets are displayed in Fig. [2] Orbital 
and physical parameters assumed in our computations are pre- 
sented in Table [2] Generally, we tried t o ma intain coherence wit h 
the data applied by iTavlor et"ai] d2007l) and iDurech et alj d2009f) . 
The effective diameter (the radius of a sphere with the same volume 
as an object) of Eger was selected indirectly: actually we scaled the 
asteroid to have the same volum e as a spheroid with semi-axes 2.3 
and 1.5 km l lBenner et al]l997l) . 

Considering the (54509) YORP, we compared two variants: a 
realistic assumption that the asteroid is an S-type object with ge- 
ometric albedo p = 0. 1, and rather ficticious case of spectral type 
E with p = 0.6. The Hapke parameters were taken from Table [T] 
If the Lambert model was assumed, the appropriate Bond albedo 



More precisely, we use the 'A-Rough' model available through the NASA 
PDS website. 
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Figure 2. Triangulated shape models: left - (54509) YORP, right - (3103) Eger. 



and emissivity from TableQ]were applied. For the Lommel-Seeliger 
model, the single scattering albedo w was computed from Ag ac- 
cording to Eq. |A3| and the emissivity followed from = 1 -A^. 
Of course, the Lommel-Seeliger scattering was not considered for 
p = 0.6 that might lead to hemispheric albedo values outside the 
[0, 1] interval. 

Figure [3] presents the simulation results for (54509) YORP. 
Although its caption mentions only Lambert and Hapke models, 
the Lommel-Seeliger results for p = 0. 1 are still there: they practi- 
cally coincide with the Lambertian solid line. We traced the values 
of (Mi) for all possible obliquities s, although the actual value for 
(5409) YORP is e = 173°. The angle between the asteroid's ver- 
nal equinox and the orbital perihelion a> , irrelevant for the rota- 
tion period related (M3}, but essential from the point of view o f 
(Mi) and (M2), responsible for the attitude tereiter et"al]|2010h . 
is o) = 102° and we used this value for all e values in Fig. [3] 
Considering (M3) values (Fig. [3] top), we observe that in spite 
of irregular shape the kind of scattering model at low p = 0. 1 has 
practically no influence on the YORP effect in rotation period, and 
even at the high albedo case (p = 0.6) the difference between Lam- 
bert and Hapke models does not exceed 10 percent. The situation 
is different for (M\) and (Mi), but there, even for the Lambert 
model, we have a dependence on albedo resulting from the heat 
conduction. Although for p = 0.1 there is almost no difference be- 
tween the Lambert, Lommel-Seeliger, and type S Hapke models, a 
high geometric albedo p = 0.6 leads to significant differences be- 
tween the Lambert approximation (dashed) and the E-type Hapke 
model (dot-dashed). The results for nominal values of e = 173° and 
ll>o = 102° are collected in Tab. [3] C omparing (M3) wit h the ob- 
served tj = (4.7 ± 0.5) x 10~ 16 rads -2 jTavlor et al .1120071) , we note 
that the present model overestim ates u> almost 7.6 ti mes, i.e. more 
than the relevant models used in iTavlor et al.1 d2007l) . On the other 
hand, the most significant part of this increase is due to the recom- 
puted reduction to the centre of mass and (more important) princi- 
pal axes system. If the original body fixed frame is used, we obtain 
a lower factor 7.0. 

In the simulations referring to (3103) Eger we compared only 
the Lambert and Hapke model for spectral type E with a high 
geometric albedo p = 0.6. In spite of a convex shape, excluding 
all shadowing effects, the dependence of all three (Mi) compo- 
nents on the scattering/emission model has the same relative mag- 
nitude as in the case of (54509) YORP. The values for the ac- 
tual spin axis orientation of Eger are provided in Tab. [4] Interest- 



Table 3. Mean YORP torques for (54509) YORP at £ = 173° and o> a = 102° . 
All values in 10~ 16 rads~ 2 . 









(Mi) 


(M 2 ) 


(M 3 ) 


Lambert 


P- 


= 0.1 


-0.21 


-0.92 


35.5 


Lommel-Seeliger 


P- 


= 0.1 


-0.18 


-0.46 


35.4 


Hapke (type S) 


P-- 


= 0.1 


-0.07 


-0.34 


34.9 


Lambert 


P- 


= 0.6 


-5.07 


-6.41 


35.5 


Hapke (type E) 


P- 


= 0.6 


-4.44 


-2.82 


32.9 



Table 4. Mean YORP torques for (3103) Eger at e = 177° and lo = 100°. 
All values in 10~ 18 rads~ 2 . 







(Mi) 


(M 2 ) 


(M 3 ) 


Lambert 


p = 0.6 


0.002 


-0.93 


1.51 


Hapke (type E) 


p = 0.6 


0.015 


-0.77 


1.39 



ingly, our modeled values of (M3) are very close to the observe d 
6j = (1.2 ± 0.8) x 10~ 18 rads" 2 reported by I6urech et all d2009h . 
Of course, this exceptional agreement can be a lucky coincidence, 
recalling inaccurate nature of photometric convex shape model, 
roughly estimated density and still a large error margin of the a> 
determination. 



5 CONCLUSIONS 

As far as photometry of Solar System bodies is concerned, the bidi- 
rectional reflectance model elaborated by Hapke leads to signifi- 
cantly different results than the basic Lambertian framework. But 
the YORP effect in rotation period occurs to be almost insensitive to 
the scattering/emission model and even at highest observed albedo 
values the difference between the two models does not increase to 
more than 10 percent. However, this low sensitivity should not be 
interpreted as the evidence of insensitivity of the scattered radia- 
tion torque on reflectivity model. Actually, the situation quite op- 
posite. Even for a given Bond albedo value, the part of the YORP 
torque originating from the recoil of reflected light significantly de- 
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Figure 3. Secular YORP effect components on (54509) YORP. Solid line - 
Lambert type S , dotted - Hapke type S (p = 0. 1 ), dashed - Lambert type 
E, dot-dashed - Hapke type E (p = 0.6). 



Figure 4. Secular YORP effect components on (3103) Eger. Solid line - 
Lambert type E , dotted - Hapke type E (p = 0.6). 



pends on the form of BRDF. But the conservation of energy implies 
that in the absence of conductivity the sum of scattered and ther- 
mally re-radiated energy is always equal to the incident energy. If 
the hemispheric albedo in some reflection model is higher than the 
Bond albedo of the Lambert case, more power is scattered, but also 
less power is thermally re-emitted and vice versa (see Fig. [5). Ac- 
tually, the same mechanism of energy balance is responsible the 
independence of YORP on albedo and emissivity in the traditional 
Rubincam approximation with Lambertian scattering/emission. In 
the ID thermal model considered in our paper the (A/3) component 
behaves exactly like in the Rubincam's approximation, so the de- 
pendence on reflectance is only due to secondary effects - mostly 
related with the small deviation of the recoil force from the normal 
to the surface. 



Using a more elaborate scattering method is more important 
in the part of the YORP effect responsible for the orientation of the 
spin axis. In Rubincam approximation, the situation is similar to 
that of (M3). But the Rubincam approximation itself is definitely 
unrealistic for the attitude even at moderate values of conductivity. 
The influence of heat conduction is proportional to the absorbed 
fraction of incident energy, hence to the albedo. It means that two 
scattering models with different dependence of hemispheric albedo 
on Sun zenith distance will differently affect the balance between 
scattered and re-radiated power. This explains why using the Hapke 
BRDF instead of Lambertian model is more important for (M\) and 
<M 2 ), then for <M 3 ). 
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Figure 5. YORP effect in a) and its components for (3103) Eger. Dashed 
line - scattered light, dotted - thermal radiation, solid line - total. Left half 
shows the Lambert case compared with the Hapke model (right half). 
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APPENDIX A: LOMMEL-SEELIGER APPROXIMATION 

The Lommel-Seeliger scattering law ( Fair bairnl2005l) is defined by 
the BRDF 



./LS(S,0) : 



(Al) 



which can be seen as a simplified Hapke model independent on 
the phase angle g. Using this simple law we find most of the ex- 
pressions in exact, closed form, depending on the single scattering 
albedo w. The hemispheric albedo is 



w / Uq 

Ah(m)=ir U+A<©ln- 

2 \ 1+// G 

leading to the Bond albedo 

2(l-ln2) 



Ab 



■w« 0.204569 w, 



and the geometric albedo 



w 



(A2) 



(A3) 



(A4) 



Note that the last formula leads to problems with w > 1 in Hapke's 
thermal radiation expressions if we try to use it for bright objects 
with p > 0.125. In these circumstances we combine the Lommel- 
Seeliger scattering with a Lambertian grey body emission model, 
imposing £ = 0. Then, 



12 



1 + 6^o + 2yu (2 + 3^ Q )ln 



A'o 



and X2 = 0. The resulting force per area 
dF 2 , vJ 

-77 =-^-[vJno + Q\n+— Xm, 

aS 3c c 



(A5) 



(A6) 



is directed along the surface normal, similarly to the Lambertian 
case. 
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APPENDIX B: DIRECT RADIATION PRESSURE 

In the usual YORP models of a single, Sun orbiting object the di- 
rect radiation pressure, opposite to the Sun vector s, is either a pri- 
ori discarded or its effect disappears after the double (rotation and 
orbit) averaging. However, this phenomenon may play some role 
when a binary system is s t udied - most notably for the BYORP 
effect dCuk & Burns! 120051 : iMcMahon & Scheeresll20ld) . In such 
cases, the force and torque should be complemented with follow- 
ing terms: 

dF d vJn Q vJjj Q 

— — = s = Qi e n + s Q m), (Bl) 

aS c c 

being the addition to Eq. J57b . and the resulting torque 

^=rA (B2) 
dS dS 

Note that for a binary system the visibility function v additionally 
involves occlusions by the second object. 

With these complements, the complete force dF c = dF + dF d 
acting on a surface element is 

dF c = - -^i[vJn Q + Q]<iS 

3 c 

+ V J. \[x i -^)dS-{X 1 +n Q s Q> )mdS], (B3) 

and the complete torque is readily obtained by the cross product 

dM c =rxdF c . (B4) 

Equation dB lb involves an implicit statement that all photons 
hitting the surface are absorbed and transfer their momentum be- 
fore being re-emitted in any form, including the one termed re- 
flection. In a perfect specular reflection, all photons arriving from 
s = //qH + s m, leave the surface in the symmetric direction s' = 
fiQit - s e m. So, the total effect of perfect specular reflection is 

dF s = —dS, (B5) 

c 

with the m component canceled. As a consequence, we may in- 
terpret the occurrence of function X2 as a footprint of imperfect 
specular reflection in the Hapke's model, with only a part of power 
incoming from s leaving the surface along s'. 



